A one-parameter family of interpolating 
kernels for Smoothed Particle Hydrodynamics 

studies 



Ruben M. Cabezon 

Dept. de Fisica i Enginyeria Nuclear, UPC. Jordi Girona, 1-3, 08034 Barcelona, 

Spain 

Domingo Garcfa-Senz 

Dept. de Fisica i Enginyeria Nuclear, UPC. Jordi Girona, 1-3, 08034 Barcelona 
and Institut d'Estudis Espacials de Catalunya. Gran Capita 2-4, Barcelona, Spain 

Antonio Relario 

Dept. de Fisica i Enginyeria Nuclear, UPC. Jordi Girona, 1-3, 08034 Barcelona, 

Spain 



Abstract 

A set of interpolating functions of the type f(v) = {sin / (f f)} n is analyzed 
in the context of the smoothed-particle hydrodynamics (SPH) technique. The be- 
haviour of these kernels for several values of the parameter n has been studied either 
analytically as well as numerically in connection with several tests carried out in 
two dimensions. The main advantage of this kernel relies in its flexibility because 
for n = 3 it is similar to the standard widely used cubic-spline, whereas for n > 3 
the interpolating function becomes more centrally condensed, being well suited to 
track discontinuities such as shock fronts and thermal waves. 
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1 Introduction 



The technique called smoothed particle hydrodynamics (SPH) was introduced 
in 1977 pQ and [2] to simulate the evolution of fluids and plasmas in three di- 
mensions. It is a gridless Lagrangian method hence one has not to be worried 
about the definition of a mesh and its further remapping because the grid is 
somehow advected by the particles themselves (see [3] for a recent review of 
this technique). A central point of the SPH formalism is the concept of inter- 
polating function (or kernel) through which the continuum properties of the 
fluid are recovered from a discrete sample of N points with mass m 8 which 
move according to the hydrodynamical laws. A good interpolating kernel must 
satisfy a few basic requirements: it must tend to a delta function in the contin- 
uum limit and has to be a continuous function with, at least, definite first and 
second derivatives. From a more practical point of view it is also advisable to 
deal with symmetric kernels of finite range, the latter to avoid N 2 calculations. 
A prototype of kernel is the Gaussian kernel: 

W G (v,h) = T7 L- 5 exp(-t; 2 ) (1) 



where d is the spatial dimension and v = |r — r'\/h is the normalized distance 
between particles in terms of the characteristic smoothing length h. The pre- 
cise value of h sets a spatial length-scale which is usually taken as the local 
resolution of SPH. Nevertheless the Gaussian kernel has an infinite range thus 
it is most practical to use a similar function but with compact support. A 
function with the required properties is the cubic M4 spline [I], [5] defined as: 
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v > 2 



(2) 



where A — |, =^ and ^ in one, two and three dimensions respectively. This 
kernel works very well and it is computationally efficient as it has been checked 
worldwide by many people [6] and [7j using SPH in many areas of physics and 
astrophysics. 

In general, hydrodynamical simulations with SPH should not significantly de- 
pend on the chosen kernel but there could be special cases in which the precise 
profile of the interpolating device becomes relevant. These situations appear 
wherever there is an abrupt change of physical variables. Trivial examples are 
strong shock formation and thermal discontinuities. In these cases the values 
of the variables around the jump region are severely damped by the interpola- 
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tion procedure and the values at the peak are underestimated. In other cases 
the numerical noise can become high enough to blur the true value of the vari- 
ables. If there are chemical or nuclear reactions, which are very sensitive to 
temperature, the outcome might be greatly altered or even completely wrong. 
Several ways to better handle discontinuities in SPH have been proposed. One 
of them is to enhance the artificial viscosity algorithm of the codes in order 
to improve the jump of the variables at the discontinuity and its thickness [8]. 
Another route is to try a more clever interpolation using kernels especially de- 
vised to handle large gradients. In particular tensorial kernels with ellipsoidal 
geometry were proposed by [9]. In normal conditions these kernels reduce to 
the spherically symmetric cubic-spline but, in the presence of a shock the 
sphere becomes an ellipsoid with its minor axis aligned with the shock direc- 
tion leading to an improvement of the resolution. However several interesting 
properties of the spherically symmetrical kernels are lost during their transfor- 
mation into ellipsoids. In particular, the use of spherically symmetrical kernels 
guarantees that any function is approximated to second order in h, thus linear 
functions are exactly reproduced. In addition energy and momentum are not 
so well conserved when ellipsoidal kernels are used. An alternative to tensorial 
interpolators is to consider strongly peaked kernels with spherical symmetry. 
In this case, for a fixed value of the smoothing parameter h, there is an in- 
crease in resolution, although that enhancement is always accompanied by an 
increase of the numerical noise [10J. In this paper we present a one-parametric 
family of spherically symmetric kernels based on harmonic-like functions (Wjf 
kernels hereafter) with compact support. Specifically we propose that the set 
of functions: 



where B n is a normalization factor, have several interesting features which 
make them suitable to SPH studies. A change in the value of the govern- 
ing parameter n lead to different shapes of the interpolative function, from 
more extended to more centrally condensed profiles as n increases. More- 
over, for n = 3 it is very similar to the well-known cubic spline given by 
Eq. (2). Thus we suggest that the use of the set of functions in current 
SPH calculations would add more flexibility to the numerical scheme with- 
out practically introducing any inconvenient. Although the function given by 
Eq. (3) is of great importance in signal analysis theory, where the case n=2 
corresponds to the so-called window function, as far as we know it has never 
been used before in connection with SPH. A somehow related interpolator 
W oc (1 — t> 2 /4)(l + cos7ru/2) was studied by [10J although only in ID. How- 
ever it was a single kernel, not a continuous family as those given by Eq.(3) 
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Interpolators of type with a high n should be compared not only to 
the M 4 but to higher order splines. In this respect the quintic M 6 polynomial 
interpolator [I] could be taken as representative. It has a compact support and 
has been used to model flows using SPH [TT]. As the M 6 spline was originally 
devised to work within a radii of 3h we have renormalized it to 2h in order to 
make plausible comparisons with W^: 



M 6 (v,h) = - 



C 
J? 



f(2 

(2- 
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(4) 
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u > 2 



where C = gnff^ anc ^ io240tt i n one ) t wo an d three dimensions respectively. 

In Section 2 we give the main mathematical features of the kernels defined by 
Eq. ([3]), and discuss their abilities to handle steep functions. A comparative 
analysis of the performance of these interpolators in disordered discrete sys- 
tems is also given in the same section. In Section 3 we analyze the behaviour 
of in connection to two classical tests: (1) the propagation of a thermal 
wave arising from a thermal discontinuity in 2D cartesian coordinates and, 
(2) the evolution of a blast wave reaching the Sedov phase, also calculated in 
2D. Finally a brief critical discussion concerning the virtues and shortcomings 
of the proposed set of functions and the main conclusions of our work is 
provided in Section 4. 



2 General properties of the set 



In this section we study the relevant mathematical features of the family of 
kernels we are proposing. Mathematical theory of interpolation tells us that 
the Fourier transform of f{w) = {sin(w)/w} k , where w = nh/ A can be used to 
generate a whole family, of polynomial interpolators of increasing degree 
[1]. For k = 2 the Fourier transform of f(w) leads to the linear interpolator 
M2 spline whereas the cubic spline is obtained for k = 4. Conversely the 
Wj? set can be regarded as the inverse Fourier transform of the My. family 
of polynomial kernels. Such reciprocal relationship basically arises because 
all acceptable kernels are merely a variation of the basic Gaussian function 
defined by Eq. ([!]). Thus the set of functions can also be considered good 
interpolating functions with the interesting peculiarity that W$ is practically 
equal to the cubic spline M 4 not after the Fourier transform of f{w) but in 
current cartesian space. 
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A plot of the profile of these kernels for several values of the governing index 
n is shown in Fig. I. As it can be seen the resulting profile for n = 3 is very 
close to that of the cubic spline, with the additional advantage that its second 
derivative is smoother. Interestingly, the profile of the quintic M 6 is well fitted 
by n = 5 (strictly, n = 4.9 provides a slightly better fit) whereas the Gaussian 
kernel is reasonably reproduced by choosing n = 2. Values of n above 4 lead to 
condensed high-peaked kernels which could be useful to handle abrupt spatial 
changes in physical variables during current hydro dynamical calculations. 



2. 1 Normalization 



Firstly we provide the normalization factor B n for in equation 3 as: 

Bn = lJ? = U (5) 



where d is the dimension of the space and, 
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(6) 



which has to be calculated numerically. In Table 1 there is shown the value of 
K n for one, two and three dimensions for different values of n. These values can 
be plotted in a diagram, Fig. 2, and approximated by a polynomial function 
which gives K n with very good accuracy for any value of n within the range 
2 < n < 7, 

K n = a 5 n 5 + a 4 n 4 + a 3 n 3 + a 2 n 2 + a x n + a (7) 



The values of the fitting coefficients Oj in ID, 2D and 3D are provided in Table 
2. The relative errors introduced by the fitting formulae are negligible, lesser 
than 10~ 5 . It is also safe to extrapolate Eq. (7) down to n = 3/2 and up to 
n = 8 without introducing a significative error. Note that the normalization 
constants K n depicted in Fig. 2 intersect around n = 7. For n > 7 the ID 
normalization constants always lie below those defined by the 2D and 3D 
lines. A plausible interpretation is that in 2D and 3D calculations with a finite 
number of particles there is a practical limit of the exponent n. Above that 
value the kernel becomes too sensible to the locus of the very firsts neighbors of 
the particle thus loosing their isotropy features. That is, the 3D kernel behaves 
pretty much as an one-dimensional interpolator. Thus we have constrained the 
value of n to be in the range 2 < n < 6 in all the numerical tests given in 
Section 3. 
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2.2 Fitting Gaussians with W 1 



Some insight about the ability of the set to reproduce steep functions 
can be gained through the study of simple Gaussians curves and surfaces 
with sharp slopes. To begin with let us consider the following one dimensional 
Gaussian density profile: 

p(v) =p (l + Re~ v2 ) ; v = x/h (8) 



which represents a Gaussian curve with a characteristic width just equal to the 
smoothing-length parameter h and maximum value p max — (1 + R)po- Thus 
the interpolators will find some difficulty to reproduce these curves because 
the resolution is similar to the width of the bell. We are mainly interested 
in three magnitudes: the maximum peak of density, the width of the curve 
at a half of that maximum and the maximum value of the density gradient. 
These Gaussians could be viewed as idealized mathematical curves mimicking 
the jump of physical magnitudes such as density or temperature associated 
to shock waves and thermal fronts. For example, taking R=3 in Eq. (jSJ) leads 
to a density jump of a factor four, the same that for a strong shock passing 
through a perfect gas with 7 = 5/3. 

From the standard SPH definitions of interpolation of a function we have (in 
ID): 



(Pma X ) = J PO (l + Re-^B. 
-2 



B n hp 





hdv 



dv + B n hp R J e v2 




dv (9) 



calling 1 1 and I2 to the integrals on the right side of Eq. ([9]) we have: 

(p max ) = Po(i + Rj^) (10) 



Therefore the closer J2//1 is to 1 the better the numerical approximation is. 
In Table 1 there is shown the fraction as a function of the parameter 

n which characterizes the harmonic-like kernels for 1, 2 and 3 dimensions. As 
it can be seen as n increases the interpolation improves substantially. For the 
particular case n = 3 the density peak does not reach the 50% of the true 
value in the 3D case whereas for n = 6 the percent rises to 64%. 
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Similarly, the maximum value of the derivative of Eq. (jSj) is: 



dp 
dx 




which takes place at v = x/h = —y2/2. The estimation of that derivative 
using the kernel is (see [3] for details on how to calculate derivatives with 
SPH): 




= np B, 

-V2/2 _ 2 
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sin 



X 



dv 



(12) 



It is easy to show that Eq. ( Tl2l) reduces to: 




-V2/2 
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(13) 



where, 




tan 



2 V 



dv 



(14) 



The factor (nl^/Ii) for the ID case as a function of n is given in the last 
column of Table 1. As we can see approaching the derivative through kernel 
estimation using Eq. (fl3i) leads to maximum values which are around the 
55% (n=3) and 71% (n=6) of the analytical estimation given by Eq. (Till . 

Finally, in Fig. 3 there is shown the complete smoothed density profile of 
the Gaussian and its spatial derivative for the cases n = 3, 5, 6 and n = 8. 
As we can see there is a clear enhancement of the resolution as n rises, in 
both the function and its derivative. However the increase in resolution is not 
monotonic because the gap between n = 3 and n = 6 is larger than that 
from n = 5 to n = 8. We can see that, in consonance to the results above, 
the width of the bell at half-height also improves as the index of the kernel 
rises. Those cases with a higher index n make also a better approach to the 
derivative of the Gaussian profile in all points, as it can also be seen in Fig. 3. 
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It is worth noting that the ability of the family to handle the above sharp 
Gaussians is independent of the particular value adopted for h (see for instance 
Eqs. 10, 11 and 13). Such good behaviour is due to the fact that both, the 
function and the kernel share the same characteristic width, h. In other cases 
the outcome of the interpolation will be dependent of the smoothing-length 
parameter. However such dependence is weak, of second order in h, [3], thus 
the discussion given in this section still holds unless the width of the Gaussian 
becomes much lesser than the smoothing-length h. 

2. 3 Pairing-instability 

Pairing instability is a well-known problem of SPH |12j . Particles that are too 
much close tend to clump together forming a stable configuration which is, 
in fact, a numerical artifact. This situation is caused by an unstable stress- 
strain relation that occurs typically when the second derivative of the kernel 
becomes negative and high strain is developed between particles. It could be 
of practical interest to have a direct control on the locus where the pairing- 
instability comes out. The higher the index n in is, the closer is that 
locus to the origin, making harder for a particle to get stuck with its closest 
neighbor. 




Solving numerically the equation (W^) (v, h) = we can find at which vq 
the second derivative becomes negative. For n = 3 we find vq = 0.6613 that 
is a bit closer to v — than the cubic spline reference value: 2/3. In Table 
3 we give the exact locus of vq for some fiducial values of n. As we can see 
the higher the exponent n is the closer vq moves to the origin, making the 
pairing-instability more unlikely to happen. 

Another interesting feature of is that the set is infinitely derivable, with 
continuous and well-behaved derivatives, while the second derivative of the 
cubic spline kernel displays several non-derivable points, as depicted in Fig. 
4. This may have relevance for those physical magnitudes that rely on sec- 
ond or even higher derivatives. Of course that shortcoming is not shared by 
higher order splines, such M$ for instance, which have well behaved second 
and higher order derivatives. In the frequent case of heat diffusion by conduc- 
tion the standard SPH implementation avoids the calculation of the second 
derivative by reducing the problem to an integral expression which uses only 
first derivatives. Nevertheless even with that formulation the smoothness of 
the second derivative may be of interest because the integral approach involves 
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a balance among the first spatial derivative calculated through the neighbours 
of a given particle. Such integral approach is somehow equivalent to perform 
second derivatives although considerably less noisy. 



2.4 Discrete disordered systems 

Interpolation in discrete systems is limited by the presence of numerical noise. 
Before attempting to use any particular kernel one has to be sure that the 
local averaged properties of the physical system are not too much distorted 
by its presence so that the noise is not taking over the dynamical evolution. 
Although to really understand the role of the disorder it would be necessary 
to solve the hydrodynamical evolution of the fluid using different levels of 
resolution, some insight can be gained by studying simple cases of interpolation 
in static lattices. Thus, to complete the study of the mathematical features of 
the proposed family of kernels, W 7 ^, we have extended the analysis given in 
the previous section to two-dimensional disordered systems. An uniform grid 
with N=57600 particles was built in a box with the masses of the particles 
conveniently crafted to reproduce the density profile of a Gaussian surface. 
As before the density jump across the bell was set to four, but this time the 
characteristic width of the Gaussian was taken large enough to ensure that 
its profile was well resolved by the SPH. The density was calculated using the 
standard SPH equation expressing mass-conservation: 



where rrij is the mass of the j-particle and n n b is the number of neighbors 
within the compact support. The number of neighbors was set constant to 
n n b = 43, meaning a ratio Ad/h = 0.305 where Ad is the typical interparticle 
separation. A level of noise was seeded by randomly re-settling these particles 
so that their final position was within a 5% radius from its initial location at 
the lattice. 

As can be seen in Fig. 5 the density profile of the Gaussian is not so much 
altered. Because the bell is wide the correct density jump, p ma x/Po = 4, is 
reached with independence of the kernel. As expected the case with n=3 led 
to a better filtering of the noise than n=5 but the differences are not large. 
Similarly there are not serious differences between the cubic-spline M4 and 
quintic Mq. Things are different when the gradient of density is computed 
using the SPH expression: 




(16) 



3=1 



(V r p); = ^TrijVrWijdTi -Tj\,h) 



(17) 



3=1 
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where V r means the radial component of the gradient. The results are sum- 
marized in Fig. 6, where it is evident that disorder has an important influence 
on the first derivative. Even though, on average, the first derivative still keeps 
with the original profile, it becomes blur showing an important dispersion, 
especially close to the peak at r ~ 10 cm. However, low-order interpolators 
M4 and W% do show a much lesser dispersion than the higher-order kernels 
M 6 and W§ , as expected. Fig. 6 also suggests that the similarity between 
the pairs M 4 , W% and M 6 , W§ kernels shown in Figs. 1 and 4 still holds in 
disordered systems, being good enough as to make them exchangeable in cal- 
culations which use a limited number of particles. A quantitative idea of the 
dispersion is given in Fig. 7 which shows the profile of the standard deviation 
a of the density gradient along the box as a function of the kernel index n. 
As we can see the dispersion becomes more important as n rises. However the 
trend is not monotonic being more accentuated for large values of n. 

It has to be stressed that in discrete systems the effect of increasing n in 
while keeping the number of neighbors, n n b, constant would somehow be 
equivalent to reduce n n b leaving the index n unaltered. For example a similar 
level of noise to that shown by n = 5 can be obtained using n = 3 but 
with n nb ~ |n n &. Nevertheless it is risky to reduce too much the number of 
neighbors because any numerical fluctuation could have a large impact on the 
smoothed variables. 



3 Numerical simulations 

The discussion given in Section 2 refers basically to the ideal mathematical 
properties of the family of interpolators. Unfortunately, as stated at the 
end of the previous section, much of these properties are partially lost in 
practical applications because of the presence of numerical noise. In current 
hydrodynamical calculations the physical system is decomposed in a finite 
number of particles with mass mi and integrals such those given by Eq. (Q or 
Eq. (TT2"j) are calculated through a summatory which involves the neighbours 
of a given particle. Therefore some level of noise is unavoidable in SPH. In 
practice a good kernel interpolator should give reliable values for the smoothed 
variables at low computational cost and keeping the noise at low enough level 
to not interfere with the simulation. It is well known that slender kernels are 
better interpolators but they also generates more noise, [TU] (see also Fig. 6). 
Thus, despite their good continuum features, choosing a too large n-value in 
Eq. ([3]) could bring more problems that solutions unless a large number of 
neighbours is used (but in that case the main advantage of increasing n is 
lost and the computational burden rises). In general the selection of the more 
suitable kernel is dependent to the particular physical situation we want to 
simulate and even to the available computational resources. 
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In this section we describe a couple of standard tests carried out in cartesian 
2D aimed at exploring the behaviour of the kernels when a limited num- 
ber of particles is used to describe the fluid. The main goal is, however, to 
provide practical examples about the advantages of using harmonic kernels 
with different index n during the simulation. The first example refers to the 
propagation of a thermal wave in an homogeneous medium. In this case the 
setting of the initial conditions has a large impact in evolution of the wave. If 
stiff initial conditions for the thermal profile are imposed we will show that the 
use of a kernel with large index n in the energy equation improves the quality 
of the simulation. The second example deals with the numerical simulation of 
a point-like explosion in an homogeneous environment, usually referred as the 
Sedov-test in the literature. In this case the use of high-peaked kernels in the 
shock front area and, conversely, low-peaked ones in the rarefaction tail leads 
to a better energy conservation. 

To perform these tests we used a 2D cartesian SPH code with temporal and 
spatially variable smoothing length. The SPH equations used in the tests cal- 
culations are the mass, momentum and energy conservation written in their 
most common formulation, [TJ. 

pi = J2 m i W v(\ r i ~ r i\' h i) ( 18 ) 



ViWuflri-rjlA^i) (19) 




dm Pi ~ 

-J7 = -^Y1 m i w H ■ vw ij(l r « - r il> K hj) 
ax Pi ■ 

"'//'./ v 'j ■ ViW u (|r< - Tj\,hi, hj) (20) 

3 



where = Vj — Vj, and the other symbols have their usual meaning. The 
kernel Wij(h{, hj) = 0.5 [Wy(/ij) + Wij(hj)} is symmetric with respect any pair 
of particles in order to ensure the exact conservation of momentum. The term 
labeled as is the artificial viscosity term defined as, 



J J IJ Vij ■ ry < 



<l,j = \ " " (2D 

Vij ■ ry > 
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and, 



P>ij 




(22) 



here £ ~ h is a characteristic length which controls the width of the shock, 
v = O.lh helps to avoid divergences when — > and the remaining symbols 
have their usual meaning. The value of the parameters a, (3 were set to 1 and 
2 respectively in all simulations. 

These equations were completed with the diffusive heat transfer equation when 
necessary: 



where q = n/(pc v ), (being k the conductivity coefficient and c v the specific 
heat), fa = 0.5(pj + pj), = — and v — O.lh is a term which avoids 
divergences when — > 0. The expression above has the peculiarity that only 
first derivatives of need to be computed at each step. 

Contributions from temporal and spatial derivatives of the smoothing-length 
parameter were neglected for simplicity. While that approach does not pose a 
problem for the first test, because the thermal wave is supposed to propagate 
through an homogeneous medium of static particles, it is not evident its valid- 
ity for the Sedov explosion calculation. It has been shown that the inclusion 
of the smoothing-length derivatives often improves the energy conservation 
|13j . although such enhancement seems only relevant for a restricted class of 
problems (for example the head on collision of two polytropes). Keeping in 
mind that the main goal of the Sedov test shown below is to make a compar- 
ative analysis among different kernels rather than to solve the problem with 
great accuracy, we have preferred not to include the derivatives of h(r, t) in 
the equations. 

Motion equations were integrated using a second order centered scheme. A 
squared lattice with one particle in each node and 1 cm of distance between 
adjacent neighbours was implemented. The box has 240 cm of side length, 
hence there are 57600 particles, and uses periodic boundary conditions. With 
only one exception (the mixed case in the Sedov test), the value of h(r, t) self- 
adapts to keep a constant number of particles, n n & = 43, within the kernel 
radius. That number of neighbors means a ratio Ad/h = 0.305 being Ad the 
interparticle separation. All particles have the same mass, adjusted to obtain 
an uniform density profile with p — 1 g.cm -3 , and obey a perfect gas EOS, 
P = (7 — l)pu with 7 = 5/3, being u the specific internal energy. 




{q i + q j ){u i -u j ){r lj ■ VWij) 



(23) 
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3.1 Thermal wave 



In this test we follow the propagation of a thermal wave moving through an ho- 
mogeneous medium with constant density. During the calculation we obliged 
the particles to be at rest so that the evolution of the wave was determined by 
evolving only the energy equation. For this particular problem the main diffi- 
culty relies in the accurate resolution of the diffusive heat transfer equation, 
Eq. (25) below. Although the numerical resolution of that equation usually 
demands the calculation of second derivatives there is a clever formulation in 
SPH, [2] that reduces the problem to first derivatives. Nevertheless such pro- 
cedure involves the balance between the first derivative of the kernel evaluated 
in discrete points inside its compact support area, which is equivalent to cal- 
culate its second derivative. As the most common used interpolator M4 does 
not have a well behaved second derivative, Fig. 4, it can be instructive to com- 
pare the evolution of the thermal wave calculated using M4 to that using the 
family and to M 6 . Another point of interest of this test relies in the hard 
initial conditions, which attempts to represent a thermal discontinuity. As in 
the case of the shock waves the multidimensional hydrocodes have also diffi- 
culties to handle thermal discontinuities which have to be artificially enlarged 
to the resolution of the code. An extreme example of a very sharp steady 
thermal wave (with a thickness between 10~ 2 — 1 cm) is the precursor thermal 
wave which induces the propagation of a self-sustained nuclear flame in Type 
la Supernovae explosions [15]. For these cases it could be very useful to use a 
highly condensed kernel, n ~ 5 — 6, in to handle the heat diffusion and 
the standard, n — 3, for the remaining equations. 

For a given initial conditions we have carried out five runs using the kernels 



3.1.1 Initial model 

We have followed the recipe by [16] (see there for details on numerical imple- 
mentation of conduction in SPH). Initializing the internal energy, following 
the two-dimensional Green's function, the system is exactly in a state that is 
a solution of the conduction problem for an initial ^-function: 



where A = 10 5 erg.cm 2 /g, a = 1 cm 2 /s is the thermal diffusivity and Uq = 10 3 
erg/g. For a given elapsed time t Eq. (124")) provides the precise profile of the 
thermal wave emerging from the initial discontinuity. At t = 0.25s the width 
of the signal becomes equal to the smoothing length h. We take the thermal 



M 4 , M 6 , Wf,Wf and W t 




(24) 
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profile at that time as the initial conditions of our simulations. That setting 
is in fact quite unfavorable for the numerical approach but it is adequate to 
show the advantages of using a kernel with a high n-value in . 



3.1.2 Results 

The evolution of a thermal wave moving inside a static media is basically deter- 
mined by the mass conservation equation, Eq. (18), and the energy equation: 

^ = - V • («VT) (25) 
at p 



where u is the specific internal energy and k is the coefficient of conductivity. 
The SPH version of Eq. (25) was Eq. (23) given above. 

The rate of change of the thermal energy content can be obtained by deriving 
Eq. (24) with respect time, 

du A ( r 2 \ 



Aat 



(26) 



The evolution of the specific internal energy profile at t=0.30 s (thus, close to 
our starting time t=0.25s), t=l s and t=5 s is shown in Fig. 8. As the initial 
state at t=0.25 s was the same for the five kernels the w-profile at the upper- 
left picture does not show significant differences. At t=l s (middle picture) the 
thermal signal has moved to the right and the peak at the origin has decreased. 
It can be seen that neither the n = 6 nor the n = 3 case correctly match the 
analytical prediction. However W§ makes a better job than W% . When the 
elapsed time was t=5 s the profile has already become smooth enough that the 
evolution was independent of the kernel index n. Nevertheless the differences 
close to the initial discontinuity (see the profiles for r < 2.5 cm) still remains. 

The differences in the outcome of the simulations can be more easily analyzed 
if we monitor the temporal evolution of the maximum value of the derivative 
of the internal energy and its position, 



du 
~dl 



(0 = T^—2 r2 ( 27 ) 



r max = 2(2at)* (28) 



The evolution of (du/dt) max is depicted in Fig. 9. Again none of the kernels 
were able to reproduce the analytical value. Such difficulty is of course caused 
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by the sharp initial profile imposed to the internal energy. Nevertheless it can 
be seen that W§ is closer to the analytical solution than W$ for the firsts 
stages of the evolution, when the width of the peak of du/dt is of the order of 
the smoothing length. The use of W§ is about a 30% more accurate than W$ 
when the system is in a such disadvantageous state for numerical simulations. 
As the system evolves the peak in du/dt widens becoming much larger than 
the smoothing length so SPH can resolve it accurately and both values of index 
n give the same results. The case with n = 5 provides a better approximation 
than that of n = 3 although slightly worse than n = 6. An inspection of Figs. 
8 and 9 also indicate that there are not significative differences between the 
evolution of the thermal wave calculated using the polynomial functions M4, 
M 6 and harmonic interpolators W%, W^ 1 . 



3.2 Sedov explosion 

In the Sedov test the evolution of a shock wave front is studied as it propagates 
in a homogeneous medium. The problem of an intense explosion in a gas is a 
standard test for hydrocodes and has some relevance in astrophysics, where is 
common to find strong shocks in many scenarios involving fluid motions at high 
velocity. The theoretical solution was found by L.I. Sedov applying self-similar 
methods and dimensional analysis for different geometries and values of 7, |17j . 
In its simplest formulation the Sedov problem has an initially cold gas at rest. 
At t=0 s there is a point explosion at the origin which in [17] was treated as an 
instantaneous release of energy at the origin and assumed that the background 
material through which the expanding gas sweeps behaves as a perfect fluid 
with P = (7 — l)pu. The most remarkable feature of this problem is that 
it leads to exact, although algebraically complicated, analytical expressions 
for the fluid variables. Unfortunately it is not easy to exactly simulate the 
evolution of the blast wave in more than one dimension. As in the thermal 
wave test the shock front was too sharp to be resolved by the hydrocode. In 
the case of SPH the artificial viscosity smears the shock over 2-3 times the 
smoothing-length. As a consequence the density jump across the shock front 
is always lesser that the factor four predicted by the theory for 7 = 5/3. 
Thus, resolution issues are here crucial not only to resolve the peak but also 
to reproduce the correct postshock variables downstream and the structure of 
the rarefied tail close to the origin. The use of adaptive kernels can greatly help 
to handle with this problem. In this respect the family of interpolators 
adds an extra degree of freedom to the scheme which, combined with a clever 
use of the adaptive smoothing-length h(r, t) is able to bring a better approach 
to the Sedov problem. 

We have carried out six calculations with the same initial conditions in order 
to analyze the influence of using kernels with different index n in the outcome 
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of the explosion. First we have simulated the evolution of the blast wave for 
W% , W§ and M 4 , M 6 comparing the results with the analytical profile. 

In our second simulation we have taken an adaptive index n in which 
changes according to the compressional state of the material. In that case a 
clear improvement of the energy conservation was seen with respect to the 
calculation that relied in the cubic spline. An alternative to Eq. (20) which 
ensures the mechanical energy conservation is to consider the evolution of 
the thermokinetic, \v 2 + u, energy [3|. Nevertheless solving the thermokinetic 
equation does not mean that the internal energy is better evolved than using 
Eq. (20), especially if there are sources or sinks of energy. 

3.2.1 Initial model 

In order to generate a shock wave, an amount of particles (about 5% of the 
total number of particles) that lie inside a box-centered Gaussian, had their 
pressure artificially raised. To smooth the initial discontinuity we take an 
initial pressure step that decays as a Gaussian function, 

2 

P(r) = P 2 + (P 1 - P 2 ) exp ^- (29) 

where Pi and P 2 are the pressures in both zones -left and right around the 
pressure step- and a sets the width of the pressure decay. In our simulations 
Pi = 10 4 dyn/cm 2 and P 2 = 1 dyn/cm 2 , to assure an internal energy reser- 
voir big enough to feed the formation of the shock wave, and a 2 = 16 cm 2 , 
which smoothes the pressure step over about two times the smoothing length. 
Again the initial profile is so sharp that SPH finds some difficulty to track 
the blast wave. In particular the conservation of energy is not longer satisfac- 
tory during the transient period until the self similar wave appears. However 
the comparison between models with different value of n is meaningful. We 
have taken the most common form of the mass conservation, momentum and 
energy equations, 

3.2.2 Results 

Six different calculations characterized by identical initial conditions but al- 
lowing variations in the smoothing kernel were carried out. Shortly after the 
induction of the initial explosion a steady self-similar state ensued. In all the 
cases we found that the evolution of the blast wave matched well the analytical 
results after the self-similar state was achieved. There were, however, inter- 
esting differences among the models which were caused by the type of kernel 
used in each calculation. The density profiles as a function of the normal- 
ized distance during the self-similar state are depicted in Fig. 10. As expected 
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M4 and W3 led to an almost identical profiles but, interestingly, the harmonic 
kernel gave a slightly better energy conservation, as it can be seen in Fig. 11. 
The simulation which used the W§ kernel led to a large enhancement in the 
energy conservation, owing to the greater weight imprinted by the nearest 
neighbors, which makes more compatible momentum and energy equations, 
Eqs. (fT9l) and fl20|) . Nevertheless the density profile during the calculation with 
n = 6 was not so well defined as in the previous cases, as it can be seen in 
the fifth snapshot of Fig. 10. There is more noise, inherited from the initial 
distribution of particles in a square lattice. Still the dispersion is affecting a re- 
duced number of particles and the profile also follows the analytical curve. The 
source of the noise can be seen in Fig. 12 which depicts several snapshots of 
the density evolution for the cases n = 3 and n = 6. In the first row, obtained 
using n = 3, the symmetry was conserved to a high degree whereas it was not 
so well preserved for n = 6 in the second row, which shows irregularities at 
0,7r/2,7r and 37r/2 rad, even at the first snapshot. Smaller irregularities are 
also apparent in the last picture at 7r/4, 37r/4, 5n/A and 7tc/A rad. These ir- 
regularities come from the small instabilities seeded by the rectangular lattice 
used to set the simulations. That geometry imposes preferred directions for 
strain propagation, a phenomenon known as hourglass instability. The more 
efficient smoothing given by low-order kernels damped the grow of that in- 
stability from the beginning in the same way that they made a more efficient 
filtering of the random noise depicted in Fig. 6. The evolution for n = 5 is close 
to that of n = 6 although a little less noisy. As in previous calculations the 
pair W§ harmonic kernel and M 6 spline led to a very similar, if not identical, 
evolution. Only the conservation of energy seems to be slightly better for W$ . 

These results suggest that the properties of low and high n-indexed kernels 
could be combined to improve the quality of the simulation without intro- 
ducing spurious numerical noise. A way to do that is to allow the smoothing 
length to take values so that the number of neighbours, n n &, of each parti- 
cle remains in a prescribed range rather than oblige them to be constant, as 
in the calculations above. Close to the shock edge the number of neighbours 
tends to rise whereas the opposite is true along the rarefaction wave. More 
(less) neighbours means less (more) numerical noise which can be compen- 
sated taking a high (low) n value in ■ Thus, in our last calculation we 
allowed n n b to stay in the range 20 < n n b < 80 and found the exponent n from 
the expression: n = 2.88539 ln(n n &) — 6.6438 which gives n=2 for n n b = 20 
and n=6 for n n & = 80. Such range in n n b led to an average in the number 
of neighbours n n b ~ 44 during the run, similar to that taken to compute the 
unmixed cases above. Knowing the particular value of n, Eq. (j7j) allows to 
compute the normalization constant. For a meaningful comparison with the 
previous simulations the parameter I used in the artificial viscosity coefficient 
given by Eq. (|22|) was kept equal to that taken to compute the n = 3, 5, 6 and 
M 4 , M 6 cases. The resulting density profile is shown in the sixth snapshot 
of Fig. 10 and the evolution in the energy conservation is depicted in Fig. 
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11. As we can see conservation of energy is half-way between that of cases 
n = 3 and n = 6 with the advantage that the numerical noise which plagued 
the n = 6 case has vanished. In addition the high-density part of the profile 
is even a little better reproduced than for the M 4 or for n = 3 case. However 
the fit of the low-density region behind the wave, although satisfactory, is not 
as good as for the n = 3 case. 



4 Discussion and Conclusions 

A one-parameter family of interpolating kernels with compact support based 
on the harmonic-like functions, W^, defined by Eq. ([3]) have been analyzed 
and checked in the context of the smoothed particle hydrodynamics method. 
Formally the widely used polynomial interpolators are linked to the proposed 
kernels via the Fourier transform of a function similar to ■ We have found 
that the set of functions defined by Eq. also displays good enough math- 
ematical features as to deserve being considered smoothing kernels by them- 
selves. Using has several interesting advantages: a) if the exponent n is 
conveniently chosen these interpolators are able to mimic with great accuracy 
various of the most common kernels used so far in SPH studies. In particular 
the cubic and quintic spline kernels M 4 and M 6 are well reproduced by taking 
n = 3 and n = 5 in Eq. (3). Even the truncated Gaussian kernel is reasonably 
fit for n = 2. Such equivalences seem to be robust in the light of the two real- 
istic test cases analyzed in Section 3, b) in the limit of a very large number of 
particles the use of functions with n > 3 improves the resolution when strong 
gradients are present (Section 2), c) Unless the cubic spline, whose second 
derivative is not smooth in several points, the has well-behaved second 
derivatives allowing the set to be derived many times, d) it adds even more 
flexibility to the SPH technique because the quality of the interpolation can 
be selected by simply varying the exponent n in equation 3. Changing the ex- 
ponent n is really straightforward and does not introduce any computational 
overload in the numerical scheme. 

In current calculations, where a moderate number of particles are put into the 
system, the choice of kernels with a large exponent n is limited by the numer- 
ical noise. Therefore the use of high-peaked interpolators should be restricted 
to special circumstances. It has been shown in Section 2 that the improvement 
in resolution is not monotonic as n rises: above n = 6 it is hard to achieve a 
much better resolution. Just the opposite is true for the noise because it grows 
faster for n > 6. Thus, it is advisable to restrict the range of the exponent 
n to the interval 2 < n < 6 in practical applications. 

Another question related to the number of particles is that either the cubic 
or the quintic spline kernels can also have their resolution increased by de- 
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creasing the smoothing-length parameter h. Reducing h has a similar effect 
as increasing n in Eq. ([3]) because there is an enhancement of the resolution 
accompanied by an increase of the numerical noise due to the reduction in 
the number of neighbours of the particle. That enhancement in the resolution 
can be used to devise an adaptive kernel scheme which has proven useful to 
handle shocks [18]. Nevertheless varying the resolution by altering h in the 
cubic spline is not as good as to change n in for two reasons. First, a 
high-peaked kernel is still able to count particles which are farther than those 
seen by a low-peaked one with comparable resolution but lesser h. In some 
circumstances the statistical weight of these particles could be high enough to 
influence the average. The second reason is simply that to control both the 
resolution and the noise it is also better to have two parameters to tune: the 
size of h, which sets the number of neighbours of a given particle, and the 
exponent n in Eq. Q. 

The only drawback of the family is that a trigonometric function has to 
be calculated every time the kernel is invoked during the simulation. That 
leads to a computational overload with respect to the evaluation of the cubic 
spline. More specifically, if we do a large number of calls to the cubic spline 
kernel and to W% with random generated arguments, the ratio in CPU time 
is roughly a factor 2.5 favourable to the polynomial kernel. Nonetheless in a 
real hydrodynamical calculation such factor is diluted by the rest of the com- 
putations especially if gravity is present or the numerical algorithm includes 
complex physics. 

As a general rule we propose to implement the kernels in the SPH equa- 
tions leaving the index n free. In normal conditions, i.e. when there are not 
strong gradients in the physical magnitudes, we should take n = 3 because for 
that value the kernel behaves as the well checked cubic spline. In several cir- 
cumstances, though, it could be wise to turn the value of the exponent n to a 
different value. It could be taken, for instance, n = 3 in all the SPH equations 
except in that one concerning the conduction transport term given by Eq. (|23|) 
if a thermal discontinuity, a hot wall for example, is present (as suggested in 
Section 3.1). When using a second order Runge-Kutta type integrators we may 
want to use a different value of n during the first and the second (centered) 
integration steps. Another possibility is to allow each particle to carry its own 
n value as in the Sedov test described in Section 3.2 where the combination 
of an adaptive n(r,t) and h(r,t) improved the fit of the density around the 
peak and led to a better energy conservation. In some cases even the choice 
of kernels with n < 3 could be worthwhile to reduce the numerical noise. 

Taking into account their easy implementation, smoothness and flexibility, 
this family of kernels are an alternative to the widely used spline kernels, 
offering several advantages that may help to improve hydrodynamical simula- 
tions which use the SPH technique. 
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Fig. 1. Profiles of W„ x h in ID for several values of n in the range n=2 to n=8. 
Superposed are also depicted the profiles of various of the most common kernels 
used in SPH: M4 (cubic-spline), Mq (quintic spline) and truncated Gaussian. 
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Fig. 3. Upper-left and right (for a zoom) : Ability of the Wjf set to fit an one-dimen- 
sional density distribution with a sharp Gaussian profile. When the characteristic 
width of the Gaussian is equal to the smoothing-length parameter h a high value of 
n leads to a better approximation to the analytical density profile. Bottom-left and 
right (for a zoom): The derivative of density is also better approached by taking a 
large n. Note that the enhancement in resolution is not monotonic as n rises. 
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Fig. 4. First (left) and second (right) derivatives of for n = 3, 5, 6 and splines 
M^jMq. The first derivative of the pair M^,W^ behaves similar but the second 
derivative profile of W% is smoother. The pair Mq,W^ shows a good match in 
both the first and the second derivatives. 
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Fig. 5. Density profile resulting from a two-dimensional distribution of slightly dis- 
ordered particles aimed at reproducing a blurred Gaussian surface. Low-order in- 
terpolators such M4 or W% (upper-left and right) give a more defined profile than 
the high-order ones Mq and W5 , but the differences are not large. 
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Fig. 6. Same as Fig. 5 but for the density gradient. Now the imprint of the disorder 
is much larger. Low-order kernels give a more efficient filtering of the noise 
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Fig. 7. Standard deviation of the density gradient profile shown in Fig. 6 for n = 2 
n = 8 and M4, Mq. The increase in the standard deviation is not monotonic 
n rises. 
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Fig. 8. Profiles of a thermal wave arising from a discontinuity for several elapsed 
times. The cubic spline and the n = 3 case give very similar results. The evolution 
calculated using the quintic spline or W§ is also very similar. The choice of a kernel 
with a high exponent n in the energy equation improves the simulation. 
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Fig. 9. Temporal evolution of the maximum flux of heat during the propagation of 
the thermal wave depicted in Fig. 8. After t=l s the evolution was independent of 
the particular index n of the kernel because the thermal profile has become quite 
smooth 
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Fig. 10. Sedov test. Density profiles during the self-similar stage for differ- 
ent choices of the index n in W^: n = 3 (upper- left), n = 5 (mid- 
dle-left) and n = 6 (bottom left). The profiles using M4 (upper-right) and 
Mq (middle-right) are also depicted. The n-mixed case is shown at the bot- 
tom right. The maximum value of the peaks are p max = 3.43 g.cm -3 (cu- 
bic-spline), p m ax = 3.42 g.cm" 3 (n = 3), p m ax = 3.53g.cm" 3 (quintic 
spline), p max = 3.52 g.cm -3 (n = 5), p max = 3.55 g.cm -3 (n = 6) and 
Pmax = 3.47 g.cm -3 (n-mixed). 31 




Fig. 11. Energy conservation during the Sedov test using different interpolators. The 
bad conservation especially shown by low-order kernels was due to the hard initial 
conditions. 
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Fig. 12. Colour-map showing the evolution of density of a point-like explosion in an 
homogeneous media (Sedov test) at times (from left to right) t=0.2 s, 0.6 s, 1.0 s 
and 1.5 s respectively. The upper-row was obtained using W$ ■ The lower row is for 
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Table 1 

Normalization coefficients K n for different values of n and dimensions (columns 2,3 
and 4). Quotient of integrals h and h in Eq. (9) (columns 5,6 and 7) related to the 
maximum value of a sharp Gaussian curve. Factor (n{h/h) in Eq. (13) (only for 
ID, column 8) related to the maximum of the derivative of the Gaussian bell. 





h/h 


n(h/h) 


n 


ID 


2D 


3D 


ID 


2D 


3D 


ID 


1 


0.424095 


0.196350 


0.098175 


0.622276 


0.417078 


0.297324 


0.263691 


2 


0.553818 


0.322194 


0.196350 


0.714339 


0.529116 


0.404909 


0.393104 


3 


0.660203 


0.450733 


0.317878 


0.769870 


0.605666 


0.486165 


0.476948 


4 


0.752215 


0.580312 


0.458918 


0.807105 


0.660911 


0.548647 


0.53.5367 


5 


0.834354 


0.710379 


0.617013 


0.833859 


0.702584 


0.597875 


0.578325 


6 


0.909205 


0.840710 


0.790450 


0.854038 


0.735124 


0.637552 


0.611219 


7 


0.978402 


0.971197 


0.977949 


0.869814 


0.761234 


0.670167 


0.637203 


8 


1.043052 


1.101785 


1.178511 


0.882493 


0.782649 


0.697430 


0.658243 


9 


1.103944 


1.232440 


1.391322 


0.892909 


0.800531 


0.720550 


0.675626 


10 


1.161662 


1.363143 


1.615708 


0.901621 


0.815690 


0.740399 


0.690227 
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Table 2 

Fitting coefficients to the normalization constant 





ID 


2D 


3D 




2.645649xl(T 


-1 


7.332473x10- 


-2 


2.719002x10- 


-2 


CL\ 


1.824975x10" 


-1 


1.196425x10" 


-1 


5.469083x10- 


-2 


a-2 


-2.426267x10' 


-2 


3.319287x10- 


-3 


1.711166x10" 


-2 


a-i 


3.112410x10" 


-3 


-5.511885x10" 


-4 


-1.237265x10" 


-3 


04 


-2.404560x10" 


-4 


4.828286x10- 


-5 


8.193975x10- 


-5 


a 5 


8.032609x10" 


-6 


-1.733766x10 


-6 


-2.552696x10" 


-6 
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Table 3 

Critical value vq at which 
the second derivative of 

and M4,Mq becomes negative 
leading to pairing-instability. 



Interpolator 


vo 


M 4 


2/3 


M 6 


0.5062 


W* 


0.6613 


w* 


0.5039 


wf 


0.4582 


w? 


0.3718 
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